{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example will setup the required electronic structures for usage in TBtrans.  \n",
    "You will also learn the importance of performing $k$-point convergence tests for systems using TBtrans.  \n",
    "We will continue with the graphene nearest neighbour tight-binding model and perform simple transport calculations using TBtrans.  \n",
    "Our example will again concentrate on graphene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(orthogonal=True)\n",
    "H = sisl.Hamiltonian(graphene)\n",
    "H.construct([[0.1, 1.43], [0., -2.7]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the above call of the graphene lattice is different from [TB 1](../TB_01/run.ipynb). In this example we will create an *orthogonal* graphene lattice, i.e. the lattice vectors are orthogonal to each other, unlike the minimal graphene lattice.  \n",
    "The minimal orthogonal graphene lattice consists of 4 Carbon atoms.\n",
    "\n",
    "Assert that we have 16 non-zero elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hamiltonian we have thus far created will be our *electrode*. Lets write it to a TBtrans readable file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H.write('ELEC.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now a file `ELEC.nc` exists in the folder and it contains all the information (and more) that TBtrans requires to construct the self-energies for the electrode.\n",
    "\n",
    "All that is required now is the device region.  \n",
    "An important aspect of ***any*** transport setup is that the electrodes *must* **not** have matrix elements crossing the device region. I.e. there must not be matrix elements between any of the electrodes. This restriction is easily accommodated in tight-binding setups, but for DFT systems it is less transparent.  \n",
    "In this tight-binding setup it simply means a repetition of the electrode 3 times; 1) left electrode, 2) scattering region, 3) right electrode.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Creating the device, `Geometry`; `Hamiltonian`; `Hamiltonian.construct`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we tile the orthogonal graphene lattice 3 times along the second lattice vector (Python is 0-based) and subsequently construct it using the same parameters.  \n",
    "This method of specifying all matrix elements is the most usable and easy scheme that is available in `sisl`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = graphene.tile(3, axis=1)\n",
    "H_device = sisl.Hamiltonian(device)\n",
    "H_device.construct([[0.1, 1.43], [0, -2.7]])\n",
    "print(H_device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Creating the device, `Hamiltonian` $\\to$ `Hamiltonian`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Geometry.tile` function is an explicit method to create bigger lattices from a smaller reference lattice. However, the `tile` routine is also available to the `Hamiltonian` object. Not only it is much easier to use, it also presents these advantages:\n",
    "\n",
    "* It guarantees that the matrix elements are the same as the reference `Hamiltonian`, i.e. you need not specify the parameters to `construct` twice,\n",
    "* It is *much* faster when creating $>500,000$ samples from smaller reference systems,\n",
    "* It also requires less code, which increases readability and is less prone to errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_device = H.tile(3, axis=1)\n",
    "print(H_device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information, you may execute the following lines to view the :\n",
    "\n",
    "    help(Geometry.tile)\n",
    "    help(Hamiltonian.tile)\n",
    "\n",
    "\n",
    "Now we have created the device electronic structure. The final step is to store it in a TBtrans readable format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H_device.write('DEVICE.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it may be convenient to plot the entries of the matrix to assert the symmetry and structure. The second line asserts that it is indeed a Hermitian matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.spy(H_device.Hk());\n",
    "print('Hermitian deviation: ',np.amax(np.abs(H.Hk() - H.Hk().T.conj())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First run of TBtrans\n",
    "\n",
    "You should first run `tbtrans` in one of the following ways (the `RUN.fdf` file is already prepared with enough input options for a successful run):\n",
    "\n",
    "    tbtrans RUN.fdf                 -->  to print output on screen\n",
    "    tbtrans RUN.fdf > TBT.out       -->  to write output on file (rename `TBT.out` if desired)\n",
    "    tbtrans RUN.fdf | tee TBT.out   -->  to print output on screen *and* write output on file\n",
    "    \n",
    "After TBtrans is complete, a number of files will be present:\n",
    "\n",
    "- `siesta.TBT.nc`  \n",
    "    The main data-file of TBtrans, this contains *all* calculated quantities, and everything that can be orbital resolved *is* orbital resolved, such as density of states.\n",
    "    \n",
    "- `siesta.TBT.CC`  \n",
    "    The energy points at which TBtrans has calculated physical quantities.\n",
    "    \n",
    "- `siesta.TBT.KP`  \n",
    "    Used $k$-points and their corresponding weights for integrating the Brillouin zone.\n",
    "    \n",
    "- `siesta.TBT.TRANS_Left-Right`  \n",
    "    The $k$-resolved transmission from `Left` to `Right` electrode. This is a consecutive list of transmissions for all energy points. Each $k$-point transmission is separated with a description of the $k$-point and its weight.\n",
    "    \n",
    "- `siesta.TBT.AVTRANS_Left-Right`  \n",
    "    The $k$-averaged transmission from `Left` to `Right` electrode. This is the $k$-averaged equivalent of `siesta.TBT.TRANS_Left-Right`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbt = sisl.get_sile('siesta.TBT.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After calculating the transport properties of the transport problem, you may also use `sisl` to interact with the TBtrans output (in the `*.TBT.nc` file):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tbt.E, tbt.transmission(), label='k-averaged');\n",
    "plt.plot(tbt.E, tbt.transmission(kavg=tbt.kindex([0, 0, 0])), label=r'$\\Gamma$'); \n",
    "plt.xlabel('Energy [eV]'); plt.ylabel('Transmission'); plt.ylim([0, None]); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several function calls present in the above code:\n",
    "\n",
    "* `get_sile`  \n",
    "is a `sisl` function to read and parse *any* file that is enabled through `sisl`. You can check the documentation to find the available files. Here we use it to make `tbt` be an object with all the information that is present in the `siesta.TBT.nc` file.\n",
    "* `tbt.transmission`  \n",
    "is a function that retrieves the transmission function from the file. It has three optional arguments, the first two being the origin electrode and the second the absorbing electrode. They are defaulting to the first and second electrode, respectively.\n",
    "`tbt.transmission` takes a third and optional argument, if `True`, or not specified, it returns the **k**-averaged transmission, else one may provide an array of integers that represent the internal k-points. I.e. the code above searches for the **k**-index of the $\\Gamma$ point, and requests only that sampled transmission function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will see a *very* crude step-like transmission function.  \n",
    "\n",
    "- Why is it not smooth, V-shaped (as it should be)? Can you change something to obtain a smooth transmission function?  \n",
    "- Why is the $\\Gamma$ transmission a fixed non-zero value? Should it be zero somewhere?\n",
    "  *HINT: check out the energies used for evaluating the transmission function*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `siesta.TBT.nc` file also contains two different density-of-states quantities. How do they differ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tbt.E, tbt.DOS(), label='DOS'); \n",
    "plt.plot(tbt.E, tbt.ADOS(), label='ADOS'); \n",
    "plt.xlabel('Energy [eV]'); plt.ylabel('DOS [1/eV]'); plt.ylim([0, None]); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sisl` is capable of interacting *much* more with TBtrans output in various ways. We will return to this in later examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Go to `In [4]` and change `H.write('ELEC.nc')` to `H.write('ELEC.TSHS')`, then adapt `RUN.fdf` file to let the left electrode read the electronic structure from the `TSHS` file instead. Redo the calculation and check if the output is the same.\n",
    "- Try and change the on-site term of _one_ of the device atoms, then redo the calculation. What is the (periodic) effect of changing the on-site for one atom? (difference should be larger than `0.5 eV`)\n",
    "- Try and create a big graphene flake with an associated Hamiltonian composed of ca. 100,000 atoms. When doing this it is essential you only `construct` the minimal graphene flake (otherwise creating such a big Hamiltonian will take several minutes). Then use `tile` or `repeat` to expand to the big system. If you try to run transport on this device, it will take forever, and may even crash. To handle such big systems, you may use the fdf-flag `TBT.Atoms.Device` to limit the calculation region (advice is to *not* do the calculation). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learned lessons\n",
    "\n",
    "- Easy creation of Hamiltonians with constant parameters (`Hamiltonian.construct`)\n",
    "- Creation of the Hamiltonian using `Hamiltonian.tile` for *very* large systems.\n",
    "- Saving the electronic structure in various formats (`Hamiltonian.write`)\n",
    "- Data extraction from `*.TBT.nc` file (e.g. `tbt.DOS`, `tbt.ADOS`, etc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
